Infection by a helminth parasite is associated with changes in DNA methylation in the house sparrow

Abstract Parasites can exert strong selective pressures on their hosts and influence the evolution of host immunity. While several studies have examined the genetic basis for parasite resistance, the role of epigenetics in the immune response to parasites is less understood. Yet, epigenetic modifications, such as changes in DNA methylation, may allow species to respond rapidly to parasite prevalence or virulence. To test the role of DNA methylation in relation to parasite infection, we examined genome‐wide DNA methylation before and during infection by a parasitic nematode, Syngamus trachea, in a natural population of house sparrows (Passer domesticus) using reduced representation bisulfite sequencing (RRBS). We found that DNA methylation levels were slightly lower in infected house sparrows, and we identified candidate genes relating to the initial immune response, activation of innate and adaptive immunity, and mucus membrane functional integrity that were differentially methylated between infected and control birds. Subsequently, we used methylation‐sensitive high‐resolution melting (MS‐HRM) analyses to verify the relationship between methylation proportion and S. trachea infection status at two candidate genes in a larger sample dataset. We found that methylation level at NR1D1, but not CLDN22, remained related to infection status and that juvenile recruitment probability was positively related to methylation level at NR1D1. This underscores the importance of performing follow‐up studies on candidate genes. Our findings demonstrate that plasticity in the immune response to parasites can be epigenetically mediated and highlight the potential for epigenetic studies in natural populations to provide further mechanistic insight into host–parasite interactions.


| INTRODUC TI ON
Parasites are important drivers of natural selection and can have major effects on individual fitness, population dynamics, and evolutionary processes in natural populations (Morgan et al., 2004).
Complex changes in parasite prevalence and virulence that may have negative consequences for species communities are expected in the current period of rapid environmental change (Altizer et al., 2013).
Resistance to parasites can have a genetic basis involving numerous genetic pathways that shape total immune response (Allen & Maizels, 2011;Hagai et al., 2018). Transcriptional regulation of immune genes has a pivotal role in mounting appropriate immune response to pathogens, through mechanisms including regulatory action of non-coding RNAs, transcriptional control of B and T Cell development (Rothenberg, 2014), and changes in chromatin structure that contribute to regulation of both innate and adaptive immunity (Smale et al., 2014). Infection by parasites is known to alter expression of immune genes, which has been well-studied in model organisms (Maizels et al., 2018) and in livestock, including sheep (Andronicos et al., 2010;Pemberton et al., 2011), cattle , and chicken (Dalgaard et al., 2015). However, the factors that modulate these changes in gene expression remain largely unknown. Thus, investigating epigenetic mechanisms that can modulate regulation of immune genes is central to understanding the molecular basis of the immune response to parasite infection.
DNA methylation is the most well-studied form of epigenetic modification, and there is mounting evidence that DNA methylation variation at loci involved in both innate and adaptive immunity are integral to the immune response to pathogens (Kondilis-Mangum & Wade, 2013;Saeed et al., 2014;Weng et al., 2012;Zhang & Cao, 2019).
DNA methylation involves the addition of a methyl (-CH3) group on the 5′ carbon of a cytosine residue by DNA methyltransferase. In vertebrates, methylation predominantly occurs at CG (CpG) dinucleotides but may also occur at non-CpG sites (Auclair & Weber, 2012).
Cytosine methylation at a promoter or transcription start site (TSS) generally reduces gene expression via inhibition of transcription factor binding or reduction of transcription rate (Deaton & Bird, 2011), whereas CpG methylation in other genomic contexts may serve a regulatory role or promote expression of non-coding RNAs (Auclair & Weber, 2012). In birds, gene body methylation does not appear to influence gene expression levels in blood (Watson et al., 2020), but may influence expression in other tissues, notably brain tissue (Derks et al., 2016).
The environmentally responsive nature of DNA methylation (Angers et al., 2010) provides a mechanism for rapid response of host's immune systems to environmental pathogens, the effects of which may be either transient or trans-generational (Poulin & Thomas, 2008;Roth et al., 2018). However, relatively few studies to date have investigated changes in DNA methylation in response to parasite infection. In mice, DNA methylation has been shown to play a central role in dendritic cell-induced Th2 immunity to helminth parasites, via regulation of key target genes (Cook et al., 2015).
Parasite-mediated selection has also been found to modulate DNA methylation in natural populations of killifish (Kryptolebias hermaphroditus; Berbel-Filho et al., 2019). In Trinidadian guppies (Poecilia reticulata) experimentally infected by an ectoparasitic monogenean, the authors identified differentially methylated regions between infected and control fish that overlapped genomic regions relevant to immune response (Hu et al., 2018). In three-spined stickleback (Gasterosteus aculeatus) experimentally infected by a parasitic nematode, genome-wide changes in DNA methylation have been observed (Sagonas et al., 2020). A study on DNA methylation in a natural population of red grouse (Lagopus lagopus scotica) also detected an epigenetic signature of parasite infection in epilocusspecific differentiation at immune genes and sites relating to histone acetylation, whereas genome-wide changes in methylation levels related to parasite load were not present (Wenzel & Piertney, 2014).
Conversely, experimental manipulation of parasite presence in mockingbirds (Mimus parvulus) did not result in a strong epigenetic signature of infection (McNew et al., 2021). Parasites may also modulate the immune response of their hosts by inducing DNA methylation changes that dampen the inflammatory response to infection (Varyani et al., 2017;Zakeri et al., 2018).
The house sparrow (Passer domesticus, Linnaeus, 1758) is a passerine bird with worldwide distribution that has been used as an ecological model species in several genetic architecture studies on morphological (Jensen et al., 2003;Lundregan et al., 2018;Silva et al., 2017), physiological (Andrew et al., , 2019, and parasite resistance traits (Lundregan et al., 2020). Individual-level data on infection by a parasitic nematode, Syngamus trachea (Montagu, 1811), have been systematically collected by feces sampling in our study system of a metapopulation of house sparrows in northern Norway (Holand et al., 2013. Syngamus trachea is circumglobally distributed and infects most terrestrial bird genera (Atkinson et al., 2008; | 3 of 34 LUNDREGAN et al. Adaptive immunity has also been shown to be important in the defense against parasitic nematodes in chicken (Andersen et al., 2013), and several poultry species have been shown to mount an adaptive immune response to infection by S. trachea (Olivier, 1944). Furthermore, S. trachea has been found to secrete proteolytic enzymes that act as immunomodulatory agents and interfere with immune signaling pathways (Riga et al., 1995), which could lead to downstream epigenetic alterations in infected house sparrows. Syngamus trachea prevalence varies spatiotemporally within our study system (Holand et al., 2013), is positively associated with temperature, and is higher following mild winters . Infection by the parasite reduces survival probability in juvenile and adult house sparrows (Holand et al., 2014) and negatively influences female reproductive success (Holand et al., 2015). Resistance to S. trachea has been shown to be polygenic in nature in the house sparrow and has been associated with several genes involved in innate and adaptive immune function, with genes linked to mucus membrane integrity and ciliogenesis, and with genes involved in physiological processes such as production of reactive oxygen species and vitamin A synthesis (Lundregan et al., 2020). The well-documented effects of S. trachea infection in this metapopulation of house sparrows, alongside genomic resources in the form of an annotated genome (Elgvin et al., 2017), make the house sparrow-nematode system a particularly well-suited study system for further research into possible epigenetic drivers of resistance to parasites in natural populations.
In this study, we use reduced representation bisulfite sequencing (RRBS) to assess whether infection by S. trachea alters genome-wide DNA methylation levels, as well as site-specific patterns of DNA methylation in house sparrows. We used a repeated sampling design to sample individuals while still in the nest (henceforth referred to as the "nestling stage") when they were assumed to be uninfected, and again after fledging (henceforth the "fledged juvenile stage") when individuals in our "case" group were infected by S. trachea. This was done to help us understand the causality of any DNA methylation differences between case and control birds: do epigenetic differences at the nestling stage when all birds are uninfected influence later probability of infection, or does infection result in subsequent DNA methylation changes?
First, we explored whether genome-wide DNA methylation profiles were different between case and control birds using redundancy analyses (RDA). Next, we used differential methylation analysis to determine whether any parasite-induced changes in DNA methylation were associated with genes relating to immune function, physical expulsion of parasites (for example by influencing mucus production, ciliary function, or mucus membrane integrity) or physiological processes relating to immunity. Then, we used gene ontology (GO) analysis to determine whether any of the CpG sites associated with parasite infection status were enriched for functional groups. Finally, we used methylation-sensitive highresolution melting (MS-HRM) analyses to validate the relationship between methylation levels at two of our candidate genes and S. trachea infection status in a larger sample dataset.

| Sampling and experimental design
The individuals used in this study were sampled as part of our longterm study of an insular house sparrow metapopulation off the coast of northern Norway (66°30′N, 12°30′E). The RRBS dataset included 12 female birds that were nestlings and fledged juveniles on Hestmannøy in 2011, and 10 female birds that were adults on one of three islands in the metapopulation (Hestmannøy, Gjerøy, or Aldra) between 2008 and 2011 (Table 1) Centre. All birds were ringed with a unique combination of three colored plastic rings as well as a metal ring with an individual identification number, either as nestlings in the nest or upon first capture by mist netting. A 25 μl blood sample was collected from the brachial vein for nestlings in the nest, and for fledged juveniles and adults upon each capture occasion. Blood samples were stored in 96% ethanol at −20°C. Blood is the most relevant tissue for studying DNA methylation changes in response to pathogens because it contains the white blood cells that orchestrate the immune response, as well as other molecules involved in immunity. However, see the "Caveats and future directions" section for a discussion of the limitations of using whole blood. Fledged juveniles and adults were sampled for feces on each capture occasion, birds were placed in a paper bag for approximately 15 min prior to blood sample collection to allow time for them to produce a feces sample. Fecal samples were collected and stored in ~1 ml of MilliQ water at 4°C until processed. Syngamus trachea fecal egg count (FEC) was quantified using the sucrose flotation method described in Holand et al. (2013). Island-specific microsatellite pedigrees are available for all birds used in this study and were used to estimate the relatedness between individuals (Araya- Ajoy et al., 2019;Billing et al., 2012).
Influence of parasite infection on DNA methylation patterns was investigated in juvenile birds using samples from 12 female house sparrows (Table 1), six of which were infected by S. trachea in their juvenile year (cases, individuals that had a FEC greater than zero when sampled at the fledged juvenile stage) and six of which were not infected as juveniles (controls, birds that were feces sampled multiple times during their juvenile year and never had a FEC greater than zero). Each bird was blood sampled twice: once at the nestling stage (10-14 days old) and later at the fledged juvenile stage when approximately 1-month-old (range: 26-37 days old, except for one individual that was sampled for the second time at 55 days old). This Note: In the juvenile dataset, cases are birds infected by Syngamus trachea (defined as birds with a fecal egg count (FEC) > 0, n = 6 individuals) on the day the fledged juvenile blood sample was collected (sample time point 2), controls are birds that were never infected by the parasite as juveniles (n = 6 individuals). In the adult dataset, cases are birds infected by S. trachea (FEC > 0) on the day of blood sampling (n = 5 individuals), controls are birds that never had FEC > 0 and were sampled for feces several times throughout their lifespan (n = 5 individuals). Lanes were spiked with PhiX to improve nucleotide diversity.

| Sequence alignment and identification of CpG sites
Prior to alignment of the RRBS sequence data, low-quality bases and adaptor contamination were trimmed using TrimGalore! 0.6.7 with default parameters. Subsequently, trimming of diversity adapters was done using the "trimRRBSdiversityAdaptCustomers.py" Python script provided by NuGEN (v 2.3). Trimmed sequencing reads were aligned against an in-silico RRBS digested, bisulfite-converted version of the P. domesticus reference genome v 1.0 (Elgvin et al., 2017) using BiSulfite Bolt aligner (Farrell et al., 2021) in single-end mode with default alignment parameters. Methylation bias plots (Hansen et al., 2012) for each sequenced individual were produced using MethylDackel v 0.3.0 to assess methylation percentage in each position of sequence reads, no trimming was performed before downstream analysis (

| Influence of parasite infection on site-specific DNA methylation
Differential methylation analyses using the RRBS data were carried out for juvenile birds using only CpG sites with at least 15% difference in mean methylation between cases and controls, resulting in 393,138 sites for the nestling dataset, and 619,626 sites for the fledged juvenile dataset. This filtering was used to prevent large deviations from a uniform p-value distribution that arise due to inclusion of many sites that are not expected to show a difference in mean methylation between study groups (see Husby, 2022). GLMMs with binomial error distribution were implemented using the R package lme4qtl (Ziyatdinov et al., 2018), to identify differentially methylated cytosines (DMCs) while accounting for relatedness between individuals and population structure (Lea et al., 2017). Pedigree relatedness between individuals was calculated using the R package NADIV (Wolak, 2012). Differential methylation between cases and controls at both the nestling and fledged juvenile stage (4530 CpG sites were analyzed for nestlings, and 8834 sites for fledged juveniles) was assessed using where y is a two-column matrix of methylated and unmethylated counts, μ is the intercept term, x G is a vector relating individuals to case-control group, β G is the group effect, and Z r is the pedigree relationship matrix. Bonferroni correction of group effect p-values with a family-wise error rate (FWER) of .05 was used to identify CpG sites that were differentially methylated between cases and controls.
To examine whether temporal change in mean methylation level between the nestling and fledged juvenile stage differed between cases and controls (8921 analyzed CpG sites), a model including interaction between case-control group and stage was fit using where y is a two-column matrix of methylated and unmethylated counts, μ is the intercept term, x G is a vector relating individuals to case-control group, β G is the group effect, x T is a vector relating samples to stage (nestling or fledged juvenile), β T is the effect of stage, Z ID is the random effect for individual repeated measurements, and Z r is the pedigree relationship matrix. Bonferroni correction of interaction effect p-values with FWER .05 was used to identify CpG sites for which there was evidence that the change in mean methylation level between the nestling and fledged juvenile stage was different between cases and controls.

| Genomic locations and functional analysis
CpG sites in the RRBS dataset were assigned to their genomic location using an annotated version of the house sparrow genome (unpublished, from Elgvin et al., 2017). The R packages "GenomicFeatures" (Lawrence et al., 2013) and "rtracklayer" (Lawrence et al., 2009) were used to extract start and end positions for exons, first introns, other introns, transcription start sites (TSS, defined as the region 300 bp upstream to 50 bp downstream of the annotated gene start), and promoters (defined as the region 3 Kbp to 301 bp upstream of the annotated gene start). BEDtools v2.29.2 (Quinlan & Hall, 2010) was then used to define overlap between analyzed CpG sites and the specified genomic features (intergenic regions, promoters, TSS, exons, first introns, or other introns) and assign sites to these genomic features (Table A3 in Appendix 1). This annotation was used to assign the DMCs that were identified using differential methylation analyses to genes for functional gene ontology (GO) analysis. If DMCs were assigned to more than one genomic feature, for example to the promoter of one gene and the TSS of another gene, both assignments were included. Only DMCs in the TSS, promoter, or first intron of genes were used for functional analyses because the effects of methylation level in these genomic features are known (Anastasiadi et al., 2018;Auclair & Weber, 2012;Deaton & Bird, 2011).

| Candidate gene verification using MS-HRM
2.6.1 | Gene regions and primer design In differential methylation analyses with limited sample size, falsepositive DMCs are likely to be common and caution should be taken when constructing a narrative around putative candidate genes (Pavlidis et al., 2012). To address this issue, we used MS-HRM analyses to verify the influence of S. trachea infection status on DNA methylation proportion at the candidate genes we identified at the fledged juvenile stage. All fledged juveniles that were infected by S. trachea at the time of sampling (those with FEC greater than 0) on five islands within the house sparrow metapopulation between 2009 and 2012 were included in the MS-HRM analyses. Uninfectedfledged juveniles (those with FEC equal to 0) that were sampled on the same islands during the same period were then randomly selected to total 350 individuals. All eight genes that had a DMC within the TSS or promoter at the fledged juvenile stage were considered for MS-HRM analyses (See Table A4 in Appendix 1). Based on our selection criteria for CpG sites included in differential methylation analyses in the RRBS datasets (sites with at least 15% methylation difference between cases and controls), a large difference in methylation was anticipated between infected and uninfected individuals. Therefore, we used bisulfite sequencing primers (BSP) that did not contain CpG sites within the primer sequence to perform the MS-HRM analyses because BSP primers are better suited to detecting large methylation differences than primers that contain CpG sites within the primer sequence (Life Technologies, 2010).
Our BSP primers were designed using Methyl Primer Express v1.0 (Thermo Fisher Scientific) and the P. domesticus reference genome (Elgvin et al., 2017 (Hamano et al., 2016), and the "Df value" was defined as the maximum absolute value of the relative fluorescence signal differences for each sample.

| Data analysis
The Df values for the methylated standard series were used to create a standard curve for each plate. These standard curves follow the non-linear regression model described in Hamano et al. (2016).
where "a" is a coefficient, M is the methylation level of each standard sample, Df is the Df value for each standard sample, and Df max is the maximum Df value for the 100% methylated standard. The "nls" function in R 4.2.1 was used to implement the regression model and calculate the estimated value of "a" as in Qi et al. (2021). Typical standard curves for NR1D1 and CLDN22 are shown in Figure A5 in Appendix 1.
Subsequently, the "predict" function in R 4.2.1 was used to generate a table of estimated methylation proportions for a range of Df values, that was then used to determine the methylation proportion of each sample. ( To investigate the relationship between methylation proportion of the amplified regions within the promoter of NR1D1 or within the TSS of CLDN22 and S. trachea infection status, we fitted linear mixed models using lme4qtl (Ziyatdinov et al., 2018) with arcsine transformed methylation proportion as the response variable ( Figure A6a,b in Appendix 1). Either infection status was included as a binary fixed factor in models, or fecal egg count (FEC) was included as a continuous covariate. Year and island were included as random factors to control for the effect of environment on methylation proportion, individual ID was included as a random factor to control for pseudoreplication because samples were run in duplicate, and the pedigree relatedness matrix was included as a random effect to control for kinship and population structure. Because rapid changes in DNA methylation levels are common in early development, particularly in altricial birds (Watson et al., 2019), we investigated whether methylation proportion at the amplified gene regions was related to age in days to evaluate the degree to which individual age could have impacted our analyses. Information on age in days was only available for the subset of individuals that were first recorded in the nest (n = 130 individuals for NR1D1 and n = 124 individuals for CLDN22).
Linear mixed models with arcsine transformed methylation proportion as the response variable, age in days as a continuous covariate, and infection status as a fixed factor were fitted for both NR1D1 and CLDN22. Year, island, individual ID, and the pedigree relatedness matrix were again included as random factors in the models. Finally, if methylation proportion at NR1D1 or CLDN22 was related to infection status or FEC, we also investigated whether juvenile recruitment was related to methylation proportion at the amplified region of that gene. We used GlmmTMB (Brooks et al., 2017) to fit a generalized linear mixed model with recruitment as a binary response variable (0 if an individual did not recruit and 1 for recruits), infection status as a fixed factor, mean methylation proportion for sample duplicates as a continuous covariate, and an interaction between infection status and methylation proportion to explore whether plasticity in the immune response generated by methylation changes due to parasitism impacted recruitment. A relationship between juvenile recruitment probability and methylation proportion was expected for infected individuals, whereas no relationship was a priori expected for uninfected individuals. Year and island were included as random factors in the model.

| RE SULTS
The mean mapping efficiency across all RRBS libraries was ~80% (Table A1 in Appendix 1), and, after filtering, 337,524 CpG sites with shared 10× coverage were available for downstream analyses in the juvenile dataset (Table A2 in Appendix 1). Mean methylation level was slightly higher in nestlings and fledged juvenile controls than in fledged juvenile cases, although this difference was not significant  Figure A1 in Appendix 1) but still explained very little of the variance in genome-wide methylation levels (RDA1 p = .015, adjusted R 2 = .005). Of the CpG sites that showed at least 15% methylation difference between cases and controls in fledged juveniles, 5455 of these sites had 10× coverage in the adult F I G U R E 1 Mean genome-wide methylation percentage in six cases and six control house sparrows that were sampled at both the nestling and fledged juvenile stage. Cases are shown in red and controls in blue. Each individual was sampled twice, once at the nestling stage prior to infection by Syngamus trachea and again at the fledged juvenile stage when case birds were infected by the parasite (n = 24 samples in total). Error bars represent standard error of the mean methylation across 337,524 shared 10× CpG sites between individuals. Mean methylation was slightly higher in nestlings and fledged juvenile controls than in fledged juvenile cases, and there was moderate evidence that individual mean methylation level decreased between the nestling and fledged juvenile stage in cases but not in controls (p-values from paired t-tests were .014 between nestling cases and fledged juvenile cases compared to .936 between nestling controls and fledged juvenile controls). dataset and could be used to determine whether the infection status of adult house sparrows was related to methylation profile at these sites (Table A2 in Appendix 1). There was, however, little evidence that adult cases had lower mean methylation levels than adult controls (unpaired t-test p-value = .558; Figure A3 in Appendix 1) and RDA ( Figure A4 in Appendix 1) showed that infection status did not contribute to the variance in methylation profile in adult house sparrows at these sites (p = .943, adjusted R 2 = −.006).

| The effect of parasite infection on methylation patterns
At the nestling stage, 56 of 4530 analyzed CpG sites showed evidence of differential methylation between cases and controls after Bonferroni correction (Figure 3a), and 120 of 8834 analyzed sites were differentially methylated between cases and controls at the fledged juvenile stage (Figure 3b). We also found support for differences in temporal change in methylation levels between cases and controls, indicated by evidence of an interaction between infection status and stage for 289 of 5410 analyzed CpG sites (Figure 3c). See Table 3 for an overview of the number of differentially methylated sites detected in each differential methylation analysis, as well as mean methylation percentage at these sites in cases and controls.

| Genomic locations and gene annotation
CpG sites shared between all individuals at the nestling and fledged juvenile stage were assigned to the promoter, TSS, exon, first intron, or other introns of genes or were categorized as intergenic.
Methylation levels were highest in exons, introns, and intergenic regions, and lowest in the TSS and promoter (Figure 4a). There was very strong evidence that methylation levels differed between all genomic features, except for introns and intergenic regions (pairwise Wilcoxon rank sum tests, see Table A3 in Appendix 1). Methylation differences at CpG sites in the TSS or promoter have been shown to influence gene expression in passerines (Laine et al., 2016), and a consistent negative relationship between DNA methylation levels in the first intron of genes and gene expression has been demonstrated in several species (Anastasiadi et al., 2018). As such, any DMCs assigned to these genomic features are of particular interest.
Of the DMCs that were differentially methylated between cases and controls at the nestling stage, 1.7% were located in the TSS, 12.8% were located in the promoter, and 0.0% were located within the first intron of genes ( Figure 4b). Of the DMCs that were differentially methylated between cases and controls at the fledged juvenile stage, 1.6% were located in the TSS, 5.0% were in promoters, and 2.5% were within the first intron ( Figure 4b). Furthermore, 3.0% of the DMCs where temporal change in methylation level F I G U R E 2 Redundancy analyses (RDA) of methylation profiles based on all CpG sites with shared 10× coverage in the juvenile dataset (337,524 sites). Each point represents one individual and each individual was sampled twice, once at the nestling stage and again at the fledged juvenile stage. RDA on DNA methylation profiles was performed separately on the nestling and fledged juvenile samples, with casecontrol group included as a fixed factor in both models. Cases are shown in red and controls in blue, and the full sibling pair is indicated by an asterisk at both stages. Case control group did not predict similarity of DNA methylation profiles between individuals at the nestling stage (p = .680, adjusted R 2 = −.001), nor at the fledged juvenile stage (p = .127, adjusted R 2 = .002), and the proportion of the variation in DNA methylation explained by case-control group at both the nestling and fledged juvenile stage was low.
was different between cases and controls were located in the TSS, 11.0% were in promoters, and 1.4% were within the first intron ( Figure 4b). Proportion of CpGs assigned to each genomic location was similar between all analyzed CpGs and DMCs. A majority of the genes that had DMCs in their TSS, promoter, or first intron have known immunomodulatory or mucus membrane integrity functions, see Table A4 in Appendix 1 for details.

| Functional analysis
Only genes with DMCs in the TSS, promoter, or first intron were included in functional analysis. No functional terms were enriched in functional analysis on the candidate genes identified for nestlings. The GO Biological Process term "epithelial cell morphogenesis" (GO:0003382) was enriched in functional analysis on the F I G U R E 3 Differential methylation analysis results. Only CpG sites that showed at least 15% methylation difference between cases (n = 6) and controls (n = 6) were analyzed, and for all models, Bonferroni correction with a FWER of 0.05 was used to identify differentially methylated cytosines (DMCs). "LGE22" is a linkage group corresponding to part of an unknown chromosome, "Z" is the sex chromosome, and "S" represents all scaffolds for which chromosome location is unknown. For (a) nestlings 56 of 4530 analyzed CpG sites were differentially methylated between cases and controls. For (b) fledged juveniles, 120 of 8834 analyzed CpG sites were differentially methylated between cases and controls. (c) Difference in temporal change in methylation level between cases and controls was also observed for 289 of 5410 analyzed CpG sites. Differentially methylated cytosines (DMCs) in the TSS, promoter, or first intron are highlighted in orange and labeled with corresponding genes.

TA B L E 3
Number of CpG sites assigned to different models in differential methylation analyses in lme4qtl, and mean methylation percentage ± SE for cases and controls separately.
genes identified for fledged juveniles when using both human and chicken gene ontologies. Three functional groups were enriched in functional analysis on the candidate genes identified in the temporal differential methylation analysis (see Table A5 in

| Candidate gene verification using MS-HRM
Using MS-HRM, we were able to investigate the relationship between methylation proportion at two of our candidate genes and parasitism by S. trachea in a larger dataset that included samples from 322 fledged juvenile house sparrows (Tables 2 and 4). There was very strong evidence that methylation proportion at the amplified region within the promoter of NR1D1 was higher in individuals infected by S. trachea  Table A3 in Appendix 1), first introns and intergenic regions (p = .290), and first introns compared to other introns (p = .037). The orange point above each box indicates the mean methylation level, and the horizontal black line represents median methylation level. (b) Genomic locations of differentially methylated CpG sites from differential methylation analyses. For all differential methylation analyses, the percentage of CpG sites assigned to each genomic feature (i.e. intergenic, promoter, transcription start site (TSS), exon, first intron, or intron) was similar between all analyzed versus differentially methylated sites.

TA B L E 4
Results of linear mixed-effects models using the MS-HRM dataset to validate the relationship between methylation proportion at the amplified regions within the promoter of NR1D1, or within the TSS of CLDN22, and parasitism by Syngamus trachea. Note: Effect size (above) and p-value (below in parenthesis) are given for the main effects of infection status or fecal egg count (FEC). Variance (above) and standard deviation (below in parenthesis) are given for the random factors pedigree relatedness (Relmat), individual ID, year, island, and the residual variance.
or FEC (β = 0.000, p = .890). However, in the subset of individuals with information on age in days at the time of sampling (Table A7 and Figure A6c,d in Appendix 1), there was moderate evidence that methylation proportion at CLDN22 was positively related to age in days (β = 0.0004, p = .023) although the effect size was small, whereas there was no evidence that methylation proportion at NR1D1 was related to age in days (β = 0.001, p = .165).
Because methylation proportion at NR1D1 was positively related to S. trachea infection status, we investigated whether juvenile recruitment probability was related to methylation proportion at NR1D1 or to infection status (Table 5)  We found that methylation proportion at NR1D1, but not at CLDN22 remained related to infection status, and that recruitment probability of fledged juveniles infected by S. trachea was positively related to methylation levels at NR1D1. This underscores the importance of performing follow-up studies on putative candidate genes, and ideally the need for functional studies on methylation levels at candidate genes (Gudmunds et al., 2022;Husby, 2020Husby, , 2022.

| Genome-wide DNA methylation patterns
DNA methylation levels were highest in exons followed by intergenic regions and introns, lower in promoters, and lowest in the TSS of genes (Figure 4a, Table A3 in Appendix 1). This is similar to methylation patterns found in previous genome-wide DNA methylation studies in passerines (Laine et al., 2016;Viitaniemi et al., 2019).
At the nestling stage, mean methylation levels were similar between cases (individuals that were later infected by S. trachea at the fledged juvenile stage) and controls (individuals that were not infected by S. trachea during their first year of life). However, some DNA methylation differences were observed at the fledged juvenile stage where infected birds had slightly lower methylation levels compared to uninfected controls, and mean methylation decreased to a greater extent in cases compared to controls between the two stages ( Figure 1). This suggests that infection of house sparrows by S. trachea could result in genome-wide DNA methylation changes.
Nonetheless, the results of redundancy analyses suggested that TA B L E 5 Results of generalized linear mixed-effects models using the MS-HRM dataset to investigate the relationship between recruitment probability and individual mean methylation proportion at the amplified region within the promoter of NR1D1.

Note:
We first ran the model on the entire dataset (n = 322 individuals) and included an interaction effect of methylation proportion and infection status to investigate whether plasticity in the immune response generated by methylation changes due to parasitism impacted recruitment. Subsequently, we ran the model separately on infected and uninfected individuals. Effect size (above) and p-value (below in parenthesis) are given for the main effects of methylation proportion and infection status on recruitment probability, as well as for the interaction effect of methylation proportion and infection status on recruitment probability. Variance (above) and standard deviation (below in parenthesis) are given for the random factor year and island.  (Holand et al., 2013), and an adaptive immune response to infection by the parasite has previously been implicated (Lundregan et al., 2020). Therefore, the epigenetic response to infection mounted by adult birds may differ from that of juveniles, and epigenetic differences due to factors such as developmental stage (Watson et al., 2019) or reproductive status at different points in the breeding season Viitaniemi et al., 2019) could also contribute to genome-wide differences in methylation in adult house sparrows.

| Genomic locations of differentially methylated cytosines
We observed differential methylation in the TSS, promoter, or first intron of immune genes at both the nestling and fledged juvenile stage and also identified several genes where temporal change in methylation level differed between cases and controls. The majority of CpG sites that were differentially methylated according to infection status were annotated to introns, exons, and intergenic regions rather than the TSS, promoter, or first intron. This result is similar to those from previous studies that performed RRBS using avian DNA from red blood cells (Pértille et al., 2017;Viitaniemi et al., 2019).
In the current study, 6.6%-14.5% of DMCs were mapped to the TSS or promoter (Figure 4b), which is similar to the proportion of DMCs mapped to these genomic features in previous studies on DNA methylation during parasite infection (McNew et al., 2021;Sagonas et al., 2020). The candidate genes that were identified in the present study had functions relating to both innate and adaptive immune response, as well as mucus membrane integrity and physical degradation of parasites (Table A4 in  Genes for which the temporal change in methylation level was different between cases and controls had a broad range of functions (see Table A4 in Appendix 1 for a complete overview). Case and control birds diverged epigenetically at genes relating to apoptosis (Nucleoside diphosphate kinase (DAMP) that is recognized by PRRs upon cell injury (Eleftheriadis et al., 2016). PRR signaling commonly activates the canonical NF-κB pathway and AKIP1 that showed evidence of temporal differential methylation regulates the rate of NF-κB nuclear translocation (King et al., 2011). Rab11A, a regulator of TLR4 transport (Husebye et al., 2010), was also identified. Temporal differential methylation analyses also identified genes related to adaptive immune processes including CD4+, T-cell, and B-cell expression and maturation (VPREB3, TPPP3, see Rosnet et al., 2004;Yang et al., 2020). These findings are in agreement with Perrigoue et al. (2008)

| Functional analysis
Results of GO functional analysis supported those from differential methylation analyses. In functional analysis of the candidate genes identified for fledged juveniles, the term "epithelial cell morphogenesis" was enriched. Many of the mesenchymal signaling pathways involved in epithelial morphogenesis may also be involved in epithelial repair after injury (Fang, 2000), and STC1 that was associated with this functional group has been related to wound closure in lung epithelium (Ito et al., 2014). In functional analysis of the candidate genes from the temporal differential methylation analysis, the functional group with the lowest p-value included terms related to TNFR1, NF-κB, and TNF signaling. These signaling pathways are interconnected and regulate immune function by inducing expression of pro-inflammatory cytokines and governing survival, activation, and differentiation of immune cells (Liu et al., 2017;Zhang & Cao, 2019). Thus, the enrichment of these functional terms suggests that processes relating to epithelial wound repair may be important in minimizing any mechanical damage of tracheal epithelium caused by S. trachea and that DNA methylation changes at immune genes may play a central role in mounting an immune response to parasite challenge.

| Candidate gene verification using MS-HRM
Subsequently, we investigated the relationship between methylation proportion within the TSS or promoter region of two of our candidate genes (NR1D1 and CLDN22) and parasitism by S. trachea using MS-HRM analyses and a dataset that included samples from 322 fledged juvenile house sparrows (Table 2). We found that methylation proportion of the amplified region within the promoter of NR1D1 was positively related to infection status (Table 4, Figure 5a), but unrelated to FEC, and was 64% greater in infected individuals (methylation proportion was 0.12 in uninfected individuals compared to 0.19 in infected individuals). The protein encoded by NR1D1 is a negative regulator of immune genes including, TLR4, IL-6, IL-1β, and Nlrp3 and is involved in NF-κB signaling and transcription of inflammation-related genes . Thus, higher methylation levels at NR1D1 in infected individuals may have occurred as part of the immune response to S. trachea to enable transcription of key immune genes. Methylation proportion at NR1D1 was unrelated to age in days in the subset of individuals for which we had age information (Table A7 in Appendix 1), so the observed relationship between methylation proportion at NR1D1 and infection status is unlikely to be an artifact of the rapid change in methylation levels that is commonly observed during early development (Watson et al., 2019). Furthermore, we found strong evidence that juvenile recruitment probability was positively related to methylation proportion at NR1D1 in infected individuals, but not in uninfected individuals (Table 5, Figure 5c,d), which is as expected because immune activation is costly (Graham et al., 2005;Råberg et al., 2009) and resource allocation to immune defense commonly trades off with other physiological processes such as growth (Tompkins et al., 2011).
Infection status did not impact juvenile recruitment probability (Table 5), which, although unexpected, is in agreement with previous work in our study system that found no relationship between house sparrow survival probability and S. trachea infection status, but instead found that survival probability was negatively related to infection severity (Holand et al., 2015). The results of the present study suggest that DNA methylation differences at relatively few immune genes could influence recruitment probability of juvenile house sparrows infected by S. trachea and that DNA methylation changes may alter the fitness costs of parasitism in our study system.
Conversely, we found no evidence that methylation proportion of the amplified region within the TSS of CLDN22 was related to either S. trachea infection status or FEC. Although we found moderate evidence that methylation proportion at CLDN22 was positively related to age in days (Table A7 and Figure A6d in Appendix 1), the effect size was small, which suggests that age differences between juvenile house sparrows were unlikely to have impacted our results.
Thus, the DMC within the TSS of CLDN22 that was identified in our differential methylation analysis on the RRBS dataset is likely to be a false-positive, despite previous studies that found that diverse pathogens (Soini, 2011) and helminth parasites (Su et al., 2011) can increase epithelial permeability by downregulating expression of Claudins. Due to constraints relating to MS-HRM primer design, we were only able to design functioning primers for two of the eight candidate genes that were identified in RRBS analyses at the fledged juvenile stage and that had DMCs within the TSS or promoter.
Nonetheless, using the results of our MS-HRM analyses we tentatively suggest that perhaps 50% of the candidate genes that were identified in differential methylation analyses on the RRBS dataset may genuinely be involved in an epigenetic response to parasitism.
This underscores the importance of exercising caution when constructing a narrative around putative candidate genes and highlights the value of performing follow-up studies on candidate genes for traits of interest.

| Caveats and future directions
Like many epigenetic studies in natural populations, we sampled blood, which comes with certain limitations (Husby, 2020(Husby, , 2022.
In our case, the use of whole blood can be considered a relevant tissue for studying methylation in relation to parasite infection because it contains white blood cells and other immune components (Scanes, 2015), and blood is useful because repeated sampling of the same individuals can easily be done. However, in birds red blood cells are nucleated (Scanes, 2015), so methylated DNA extracted from whole blood originates from a mixture of red and white blood cells alongside cell-free DNA and other important immune components, as well as other cell types that can interfere with methylation estimates. As there is heterogeneity in the DNA methylation profiles of different blood cell types (Adalsteinsson et al., 2012), it is possible that the changes in methylation patterns that we observed in infected juvenile house sparrows are due to changes in the cell composition of whole blood during infection, rather than due to changes in DNA methylation levels in response to immune challenge per se (Husby, 2020(Husby, , 2022. One solution to disentangle this is to examine methylation patterns in white blood cells only in future studies.
However, with present methods large blood samples are needed to separate the buffy coat that contains the white blood cells and accounts for <1% of a whole blood sample, so repeated white blood cell sampling of small passerines is not currently feasible.
Because RRBS is a reduced representation technique, this study characterized only a limited number of the approximately 15 million CpG sites in a bird genome (Derks et al., 2016). Although RRBS data are enriched for CpG-dense regulatory regions of the genome where DNA methylation changes are more likely to influence gene expression (Gu et al., 2011), it is possible that additional effects of infection by S. trachea on DNA methylation occurred in areas of the house sparrow genome that were not sequenced. Prior to differential methylation analyses, we removed CpG sites with less than 15% difference in mean methylation level between case and control F I G U R E 5 Results of MS-HRM analyses. (a) Mean methylation proportion ± SE according to infection status (uninfected = blue, infected = red) at the amplified region within the promoter of NR1D1. (b) Mean methylation proportion ± SE according to infection status (uninfected = blue, infected = red) at the amplified region within the TSS of CLDN22. (c) There was strong evidence that juvenile recruitment probability was positively related to methylation proportion at NR1D1 in infected individuals (β = 2.097, p = .007). (d) There was little evidence that juvenile recruitment probability was related to methylation proportion at NR1D1 in uninfected individuals (β = 2.374, p = .103).
birds because filtering out such sites is statistically convenient and relatively common (see e.g. Hu et al., 2018;Metzger & Schulte, 2018). However, the consequences of such filtering are not well understood and may lead to p-value inflation and overestimation of the relative importance of DNA methylation on the phenotype (Husby, 2022). Thus, future developments of methods to control for the highly skewed p-value distributions in differential methylation analyses without a priori filtering would hopefully mitigate this problem (see Husby, 2022 for more information). We were also limited by relatively low sample size in our RRBS analyses ( (Lundregan et al., preprint), and that changes in DNA methylation can lead to changes in RNA expression .

| Conclusions
Parasites are major drivers of ecological and evolutionary processes in natural populations, and exert strong selective pressures on their hosts (Altizer et al., 2013;Morgan et al., 2004). There is increasing evidence that epigenetic modifications play an important role in mounting the immune response to parasite challenge. However, many previous studies used animal experiments (Cook et al., 2015;Hu et al., 2018;Sagonas et al., 2020), Genes that were differentially methylated between cases and controls at the nestling stage were related to immune homeostasis and initial immune response, which may suggest that regulatory differences in these processes could make some birds more susceptible to parasite infection. Several genes that were differentially methylated between infected birds and controls at the fledged juvenile stage, as well as genes identified in the temporal differential methylation analysis, were related to innate and adaptive immune processes.
Thus, parasite infection may result in DNA methylation changes at diverse immune genes. Nonetheless, in differential methylation analyses with limited sample size, and in genome scan studies more generally, caution should be exercised when constructing a narrative around identified genes (Pavlidis et al., 2012). Thus, genes close to the DMCs identified in the present study should be regarded as candidate genes until verified. As a first step toward this, we used MS-HRM analyses and a larger dataset that included 322 fledged juvenile house sparrows from five islands in the metapopulation to verify the relationship between methylation proportion at NR1D1 and CLDN22 and S. trachea infection status. We found that methylation proportion at NR1D1, but not at CLDN22, remained related to infection status, which underscores the importance of performing follow-up studies on candidate genes. The observed positive relationship between juvenile recruitment probability and methylation proportion at NR1D1 in infected individuals suggests that birds may mount an epigenetic response to parasitism, which can result in a fitness advantage through increased survival probability. Taken together, the results of the present study highlight the potential for ecological epigenetics studies to provide a mechanistic understanding of host-parasite interactions in natural populations.

ACK N OWLED G M ENTS
We thank students and fieldworkers for their data collection efforts and Grethe Stavik Eggen and Randi Røsbak who provided valuable help with analyzing feces samples. We are also grateful to the inhabitants of our study area at Helgeland, whose hospitality made this work possible. We thank two anonymous reviewers for helpful feed-

CO N FLI C T O F I NTE R E S T
The authors have no affiliations or associations, professional, financial, or otherwise, that could influence our objectivity.

DATA AVA I L A B I L I T Y S TAT E M E N T
The phenotypic and epigenetic data used in this study is available on Promoter efna5b: Ephrin-A5b Region-specific apoptosis during early brain development Park et al. (2013) Note: Site position, methylation percentage in cases vs. controls, p-value, genomic region, gene, and any known gene functions are given.

TA B L E A 4 (Continued)
TA B L E A 7 Results of generalized linear mixed-effects models using the methylation-sensitive high-resolution melting (MS-HRM) dataset to investigate the relationship between methylation proportion at NR1D1 and CLDN22 and age in days Note: Models were run using the subset of individuals for which information on age in days was available (n = 130 individuals for NR1D1 and n = 124 individuals for CLDN22). Effect size (above) and p-value (below in parenthesis) are given for the main effects of infection status and age in days on methylation proportion. Variance (above) and standard deviation (below in parenthesis) are given for the random factors pedigree relatedness (Relmat), individual ID, year, island, and the residual variance.

F I G U R E A 1
Redundancy analyses (RDA) of methylation profiles based on all CpG sites with shared 10× coverage in the juvenile dataset (337,524 sites). Each point represents one individual: Cases are shown in red and controls in blue, and each individual was sampled twice, once at the nestling stage (points) and again at the fledged juvenile stage (triangles). An RDA on methylation profiles was performed on all samples, with case-control group and stage (nestling or fledged juvenile) included as a fixed factors in the model. (a) Case-control group contributed more to variation in RDA1 than stage, whereas stage contributed more to variation in RDA2. There was evidence that RDA1 was related to genome-wide methylation levels (RDA1 p = .015, RDA2 p = .562). However, due to pseudoreplication (the same individuals were sampled at both the nestling and fledged juvenile timepoints) the significance of this result should be interpreted with caution, and the proportion of the variation in methylation levels explained by case-control group and stage was low (adjusted R 2 = .005). (b) the first constrained axis (RDA1) plotted against the first unconstrained axis (PC1). Here, samples cluster by case-control group as well as by stage, although PC1 contributed more to the variance in genome-wide methylation levels than RDA1.
F I G U R E A 2 p-Value distributions from differential methylation analyses in lme4qtl. For differential methylation analysis we selected sites with at least 15% difference between cases and controls for the nestling and fledged juvenile timepoints, or with at least 15% Δ methylation difference between cases and controls in the temporal analysis.

F I G U R E A 3 Mean methylation percentage in adult samples at
CpG sites that showed at least 15% methylation difference between cases and controls in fledged juveniles. Cases are shown in red and controls in blue. Error bars represent standard error of the mean methylation between individuals across 5455 sites that had shared 10× coverage between all adult samples. There was little evidence that mean methylation across these sites was lower in adult cases than adult controls (p-value from an unpaired t-test = .558).

F I G U R E A 4
Redundancy analyses (RDA) of methylation profiles based on the 5444 CpG sites that showed at least 15% methylation difference between cases and controls in fledged juveniles and that had at least 10× coverage in the adult dataset. Each colored point represents one adult individual, infected individuals are shown in red and controls in blue. Infection status did not predict similarity of methylation profiles at these sites in adult birds (p = .943) and explained little of the variation in methylation levels between individuals (adjusted R 2 = −.006).

F I G U R E A 5
Typical standard curves for (a) NR1D1 (nuclear receptor subfamily 1 group D member 1) and (b) CLDN22 (Claudin 22). Df is the maximum absolute value of the relative fluorescence signal differences for each sample, based on the raw fluorescence data from methylation-sensitive high-resolution melting (MS-HRM).